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We develop a general formalism to treat, in general relativity, the nonradial oscillations of a 
superfluid neutron star about static (non-rotating) configurations. The matter content of these stars 
can, as a first approximation, be described by a two-fluid model: one fluid is the neutron superfluid, 
which is believed to exist in the core and inner crust of mature neutron stars; the other fluid is a 
conglomerate of all charged constituents (crust nuclei, protons, electrons, etcetera). We use a system 
of equations that governs the perturbations both of the metric and of the matter variables, whatever 
the equation of state for the two fluids. The entrainment effect is explicitly included. We also take 
£N| ' the first step towards allowing for the superfluid to be confined to a part of the star by allowing for 

an outer envelope composed of ordinary fluid. We derive and implement the junction conditions for 
the metric and matter variables at the core/envelope interface, and briefly discuss the nature of the 
involved phase-transition. We then determine the frequencies and gravitational- wave damping times 
for a simple model equation of state, incorporating entrainment through an approximation scheme 
which extends present Newtonian results to the general relativistic regime. We investigate how the 
quasinormal modes of a superfluid star are affected by changes in the entrainment parameter, and 
unveil a series of avoided crossings between the various modes. We provide a proof that, unless 
the equation of state is very special, all modes of a two-fluid star must radiate gravitationally. We 
also discuss the future detectability of pulsations in a superfluid star and argue that it may be 
possible (given advances in the relevant technology) to use gravitational-wave data to constrain the 
parameters of superfluid neutron stars. 
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pressure p-modes, and perhaps also the gravity g-modes and the Coriolis restored r-modes, are excited. By combining 
observational data with theoretical models researchers have been able to infer details of the Sun's internal structure, 
e.g. the sound speed at different depths. These impressive results of so-called "helioseismology" provide inspiration 
for further research into stellar oscillations and the hope that "asteroseismology" j| will help further unveil details 
of the structure of distant stars. The purpose of this paper is to advance our modelling of non-radial oscillations of 
non-rotating, old and cold neutron stars that contain superfluid components. 

Within the framework of general relativity oscillating compact stars provide an interesting potential source of 
gravitational radiation. There are several scenarios in which one would expect a neutron star to pulsate wildly, 
e.g. following its formation in a gravitational collapse, and it is interesting to ask whether the associated gravitational 
waves could be detected on Earth. Should this be the case, one can hope to use the gravitational-wave data to probe 
the star's interior and possibly put constraints on the supranuclear equation of state j|, [|. 

In the last few years much attention has been focussed on gravitational-wave driven mode-instabilities in rapidly 
rotating stars. Ever since the seminal work of Chandrasekhar, Friedman and Schutz |^|, fj], [| in the 1970s it has been 
known that gravitational waves can drive various modes of oscillation unstable, and recently it has been shown that 
the r-modes are particularly susceptible to the instability || |10| . Most astrophysical r-mode scenarios regard hot 
young neutron stars, but it has also been proposed fl~i| that the instability may operate in mature accreting neutron 
stars, e.g. in Low-Mass X-ray binaries. These stars are expected to have core temperatures well below the superfluid 
transition temperature, and hence one would need to account for superfluidity in any detailed r-mode model. This is 
particularly important since mutual friction may provide a strong dissipative mechanism on any mode in a superfluid 
star [ p^| . In addition to this, it has been argued that the presence of a viscous boundary layer at the base of the 
crust of a mature neutron star will lead to a very strong damping on the r-modes p3| . While the issue of mutual 
friction has been discussed by Lindblom and Mendcll [FTI) , there are as yet no studies of dissipation due to a "realistic" 



Ever since the realisation that Cepheids are stars undergoing radial pulsation , the oscillation of stars has been 
an important research area. With much improved sensitivity, observations over the last few years have established a 
plethora of pulsating stars. The best case by far is provided by the Sun for which it is known 0] that many high order 
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core-crust interface, that is likely to involve the transition from a regime where superfluid neutrons co-exist with a 
lattice of crust nuclei to a region composed of superfluid neutrons and (possibly superconducting) protons. 

Further motivation for our current work is provided by results from attempts to model the equation of state at 
supranuclear densities. Many modern equation of state calculations (within relativistic mean field theory) predict 
a sizeable hyperon fraction at high densities. However, the hyperons act as an effective "refrigerant" which means 
that the star would cool very fast. In fact, it has been argued that the presence of hyperons would lead to core 
temperatures well below those indicated by observations. The favoured resolution to this problem corresponds to the 
hyperons being superfluid. In addition, neutron stars may have cores composed of deconfined quarks which may also 
pair into exotic superfluid states. In other words, if we want to understand the dynamics of the core of a realistic 
mature neutron star we have to allow for one (or more) partially decoupled superfluid components. 

The need for an improved understanding of these problems was recently emphasized by the suggestion that the 
presence of hyperons in the deep core would lead to a strong bulk viscosity wich could potentially completely suppress 
the unstable r-modes This suggestion brings may difficult issues to the top of the agenda. For example, 

a detailed study of unstable modes must necessarily account for the fact that a superfluid component (s) may move 
relative to the normal fluid component. One might expect this to have a significant effect on the estimates of 
time scales for an unstable mode. Furthermore, as has been pointed out by Andersson and Comer ]l7[ , one would 
expect new classes of r-modes (and indeed other inertial modes) to exist in a superfluid star. This means that the 
superfluid problem is likely to be richer than the ordinary, single fluid, one. Given that there have not yet been any 
detailed investigations into these issues it may be premature to draw any definite conclusions regarding the effect that 
superfluidity may have on mode-instabilities. 

It is clear that a considerable amount of work on stars with one, or several, superfluid components remains to be 
done before we can claim to have an understanding of possible instabilities in mature neutron stars. These problems 
provide strong motivation for our current work. 

Studies of oscillating superfluid stars were pioneered by Epstein ]T3], who was the first to suggest that there ought 
to exist modes of oscillation that are unique to the superfluid case. These modes have since been calculated both 
in Newtonian theory |2(], |2l| and general relativity [^2| f23j| . It is now well established that a simple two-fluid 
model of a superfluid neutron star core has two families of fluid pulsation modes. The first of these is essentially the 
standard pressure p-modes, for which the two fluids tend to move together. The superfluid modes, on the other hand, 
are distinguished by the fact that the two fluids are counter- moving 17, 22] . Another distinguishing characteristic of 



the superfluid modes is that their mode frequencies have a fundamental dependence on entrainment. Andersson and 
Comer |17| have used a local analysis of the Newtonian superfluid equations to show that the modes are essentially 
acoustic in nature, with frequencies that depend on the stellar parameters as 

^ ^/(Z + l) 
m* r A ' 

where c p is (roughly) the speed of sound in the proton fluid, m p and m* are the bare and effective proton masses, 
respectively, r is the radial coordinate, and I is the index of the relevant spherical harmonic Yi m (6, </>). These results 
have recently been confirmed by detailed mode-calculations |2l|l . The strong dependence of the superfluid modes on 
parameters that are difficult to constrain otherwise suggest a strategy that may be used in conjunction with future 
observations to narrow down the theoretical uncertainties f23| . Specifically, an observational determination of the 
mode frequencies of an oscillating neutron star could be used to constrain the proton effective mass in dense nuclear 
matter. This would have immediate implications for BCS energy gap calculations since the effective mass is part of 
the required information (see, for instance, Khodel, Khodel, and Clark Q). 

In this paper we concentrate on the equations that describe polar perturbations (even parity) of a general relativistic 
superfluid neutron star, using as our starting point the formalism developed by Comer, Langlois, and Lin |22j . The 
neutron star is assumed to consist of two distinct regions: a core consisting of matter in superfluid states, and an 
envelope of ordinary fluid matter. We report progress in three important directions. Firstly, we have obtained the first 
ever results for gravitational-wave damping rates of the superfluid oscillation modes. This is important as it allows 
us to assess the relevance of these modes for gravitational- wave astronomy. Secondly, we have developed a framework 
in which a superfluid region can be matched to an ordinary fluid region. (Lindblom & Mendell p9[ have employed 
a similar model in the Newtonian regime.) The introduction of an ordinary fluid envelope is the first step towards 
allowing for the presence of superfluids that are confined to distinct regions in the star. However, it is important to 
point out that since we do not incorporate elasticity in our model the envelope does not play the same role as the crust 
of a mature neutron star. In essence, our model corresponds to a star in which the fluid degrees of freedom contain 
a conglomerate "proton" fluid (consisting of nuclei and electrons in the envelope, and electrons and superconducting 
protons in the core) that extends from the surface to the center of the star and superfluid neutrons that exist only 
in the core. We work out the relevant junction/boundary conditions at the core/envelope interface and determine 
the quasinormal modes for such a neutron star model. Thirdly, we have included a simple model for the so-called 
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entrainment effect, based on an expficit equation of state and an approximation that assumes that the fluid velocities 
are small compared to the speed of light. This allows us to provide the first detailed results describing the effect that 
entrainment has on the oscillation modes of a superfluid star. As we will show, this leads to the presence of so-called 
"avoided crossings" between neighbouring modes as the strength of entrainment is varied. 

The main body of the paper begins with a discussion, in Sec. [0], of the basic formalism used to describe non-radial, 
linearized oscillations of general relativistic superfluid neutron stars. This contains a discussion of the background, 
spherically symmetric and static equilibrium configurations, and the introduction of the variables that are used to 
model the non-radial oscillations. We also describe the key details of the numerical methods used to solve the linearized 
Einstein/superfluid field equations. The equation of state and inclusion of the entrainment effect are described in 
Sec. and the results of our numerical analysis are presented in Sec. [V. In Sec. |y| we prove that all pulsation 
modes of a superfluid star must radiate gravitational waves, unless the equation of state belongs to a very special 
class. Finally, we consider in Sec. [vj the question of whether or not superfluid modes, excited, for instance, during a 
puls ar gl itch, can be reliably extracted from gravitational-wave data. Some final remarks are offered in the concluding 
Sec. VII. For clarity of presentation and emphasis of the main physical results, we have relegated some of the formal, 
mathematical details to appendices. Appendix I is devoted to a derivation of the junction conditions used to tie 
together the core with the envelope. In Appendix II we present an analytical equation of state that can be used to 
incorporate the lowest-order effects of entrainment. Finally, in Appendix III a new method is presented which can be 
used to accurately determine long-lived quasinormal modes. We will be using geometrized units (except in Appendix 
II where the speed of light is restored) and "MTW" j25l conventions throughout the paper. 



II. PERTURBATIONS OF SUPERFLUID STARS 



In this Section we summarise the equations that govern the equilibrium configurations of the superfluid core and the 
normal fluid envelope. What must be determined in the core are the neutron and proton number densities, and two 
metric coefficients. In the envelope only the proton number density remains but there are still two metric coefficients 
to specify. We also describe the complicated set of coupled linear perturbation equations that need to be integrated, 
and outline the computational strategy that we have adopted. 

In what follows, the center of the star is at radial coordinate r = 0, the core/envelope interface will be at r = R c , 
and the surface at r = R. 



A. General Relativistic Superfluid Formalism 

The formalism to be used here, and motivation for it, has been described in great detail elsewhere [^2], |2^. ||^, 
[^,|30[|3l], |2||3§. Hence, we will only mention the key features here. The central quantity of the superfluid formalism 



is the so-called master function A. It depends on the three scalars n 2 — —n^n 11 , p 2 — —p^p^ and x 2 — —p^n^ that 
can be formed from the conserved neutron and proton (p^) number density currents. The master function is 
such that, when the two fluids are co-moving, — Mn 2 , p 2 , x 2 ) corresponds to the total thermodynamic energy density. 
Once the master function is provided (see Sec. |V| for a simple analytic equation of state) the stress-energy tensor is 
given by 

^ = $^+/ X , + ^, (2) 

where 

* = A-nVp-p p X P (3) 

is the generalised pressure, and 

ju„ = Bn u + Ap v , (4) 
Xv = An v + Cp u , (5) 

are the chemical potential covectors. We have also introduced 

■4 = -!^ B=-2£, C = - 2 |4. (6) 
ox z on z op 1 

The momentum covectors /i„ and %„ are dynamically, and thermodynamically, conjugate to n v and p v . The two 
covectors also make manifest the so-called entrainment effect, that affects the dynamics of a superfluid neutron star 
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in a crucial way. It is easy to see that the momentum of one constituent (p, v , say) carries along some of the mass 
current of the other constituent if A ^= (since n u is a linear combination of n v and p v ). On the other hand, if A = 0, 
i.e. if the master function does not depend on x 2 , then there is no entrainment. 

While the standard one-fluid problem can be expressed solely in terms of the Einstein equations, with the fluid 
equations of motion being automatically satisfied "by virtue of the Bianchi identities" , the two-fluid problem is 
different. Because of the additional dynamic degrees of freedom associated with the second fluid we need to use also 
(a subset of) the fluid equations of motion. The equations that need to be solved, in addition to the Einstein field 
equations, are two "continuity" equations 

VX = = . (7) 

These equations represent conservation of the superfluid neutrons and the protons separately. That is, we ignore 
"transfusion" from one component to the other due to, for example, weak interactions j3lj ]. This is likely to be a 
reasonable approximation for the timescales and amplitudes that are relevant for nonradial pulsation. We also have 
two Euler-type equations 

r^V^] = ifVbXv] = • (8) 



B. Equilibrium models 

In order to determine the background fluid configuration we need to evaluate the associated metric. We take our 
equilibrium configurations to be spherically symmetric and static, so the metric can be written in the Schwarzschild 
form 

ds 2 = -e v dt 2 + e x dr 2 + r 2 (d9 2 + sm 2 6d<j) 2 ) . (9) 

The two metric coefficients are determined from two components of the Einstein equations, which can be written in 
the form 

1 — e A 1 — e x 

A' = 8?rre A A , v 1 = h 8vrre A * , (10) 

r r 

where a primes represents a radial derivative, and it is to be understood that A = K(n,p) and ^ = ty(n,p) in the 
interval < r < R c and A = A(p) and ^ = *(p) in the interval R c < r < R. 

The e qua tions that determine the radial profiles of n(r) and p(r) in the core have been derived in Comer, Langlois, 
and Lin p^| . They follow from (|J), and can be written 

Alp' + By + hj5n + A P y = , C$p' + A%n' + ^(Cp + Any = , (11) 



where 



, n dB n dA 2 n dA 2 OA 

ap z on z ap z ox* 



»o „ n dB 2 A dA dA 2 
B° = B + 2—n 2 + A— I np+— I p 2 , 
on* on ox* 

In the coefficients above, one sets x 2 = np after the partial derivatives are taken for the equilibrium configuration. 

For the particular model equation of state we use in our numerical calculations, the above equations require that 
the two fluids have a common surface (i.e. n — > as p — > 0) if one assumes chemical equilibrium. The oscillations of 
such models were studied in p^ ] . Here we want to consider models such that the two-fluid regime is enclosed within a 
single fluid envelope. It is easy to build such a model by relaxing the assumption of chemical equilibrium in the core. 
The resultant models are likely not too unrealistic given the relatively long timescales required for nuclear interactions 
to reinstate equilibrium. 

Our models are such that the superfluid neutron number density is zero in the envelope, and thus the only matter 
equation is 

C V + \cpv' - , (13) 
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where now 



(14) 



There are three sets of "boundary" conditions that must be dealt with: a set at the center, one at the interface, 
and the remaining one at the surface of the star. In view of Eq. ([l^), requiring a non-singular behavior at the center 
of the star will impose that A(0) = 0, and consequently A'(0) and z/(0) must also vanish. This in turn implies, in 
view of Eq. (|ll|), that p'(0) and n'(0) have to vanish as well. Although our analysis of the junction conditions at the 
interface allows for other possibilities, our numerical calculations assume that the phase transition that leads to the 
formation of superfluid neutrons is second order. In other words, we consider the energy density to be continuous at 
R c . From the analysis presented in Appendix I, we see that the two metric coefficients and the pressure must then also 
be continuous. We will also assume that the proton number density and the proton chemical potential are continuous 
at the interface. At the surface, we will only consider configurations that satisfy p{R) — 0. A smooth joining of the 
interior spacetime to a Schwarzschild vacuum exterior at the surface of the star implies that the total mass M of the 
system is given by 



M = 



-4n / r 2 A(r)dr 
Jo 



(15) 



and that = 0. 



C. The linearized field equations 

It is well-known that all non-trivial pulsation modes of a nonrotating fluid star correspond to polar pertubations 
(often referred to as "even parity"). In the so-called Regge- Wheeler gauge Jm), the corresponding metric components 
are 



e"r l H {r) iur l+x H x (r) 

iwr l+1 Hx(r) eViJ (r) 

r l+2 K(r) 

r l+2 sm 2 6K(r) 



Pi{0) 



(16) 



where Pi (9) are the Legendre polynomials. This decomposition will be applied to both the core and the envelope. 

Writing the neutron and proton four-currents as — nu^ and p^ — pv^, where u^u^ = — 1 and v^v 11 = — 1, one 
can show that in the core the independent components are 



Su^e-^pC , 6v l = e-^ 2 ^Se p , 



dt 



(17) 



where 



SC n = e-^ 2 r l - l W n {r)Pie^ , 5& = - 1 *-*V n (r)-^P l e'** , 



(18) 



and 



SC P = e- x / 2 r l - 1 W p {r)P l e 1 ' 



5e p = -r l - 2 V p (r)-^Pi< 



(19) 



The Lagrangian variations in the neutron and proton number densities can be written as, respectively, 

— = Sn + n'e-^r'^Wn , 
n 

^ = Sp + p'e-^r^Wp , 
P 

and the conservation equations (Q) for the particle number currents yield 



An 
n 



-r l e- x ' 2 



l + l 



-WL 



1(1 + 1) 



1 



(20) 



(21) 
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Ap 

P 



= -r l e- x ' 2 



l -^w P + \w' p 



Kl + l) y 1„ K 



(22) 



Thus, all matter variables can be expressed in terms of the velocity variables W n , p and V n ,p- In the envelope, we have 
as the only independent matter variables W p and V p . 

The set of perturbation equations that we solve for in the superfiuid core have already been listed in j2^] , but since 
our core/envelope problem requires a slightly different method of solution we repeat the relevant equations here. We 
also need these equations for the analysis in Section V. 

First we define (in analogy with Lindblom and Detweiler's |35[ ^6) approach to the one-fluid problem) two new 
variables as 



and 



X„ = 



X P =p 



W2 



liH + e""/ V (BnV n + ApV p ) 



l"-M!L(B°nW n + J%pW p ) 



v/2 



X H a + e- u ' 2 u 2 (CpV p + AnV n ) 



>-A)/2P_ (M) 



{C ( > P W P + A° nW n 



Then we find that the Einstein and superfiuid field equations yield an algebraic constraint equation 



2-1 -I 2 



*(l- e -A)-ftr* 



-2e x ~ v u? + e x 



I 2 +1-2 



-,2A 



H + 

1 - e~ x 



2uj 2 1(1 + 1) x fl-e 
~ 2 6 



— A 



8tt* 



+ 8tt* ) ( 1 (l - e~ A ) - 4irr 2 ty 



Hi 
K 



+\^e x - vj2 {X n + X p ) = , 
and a system of coupled ordinary differential equations (where we use the definition T>q = BqCq — (_4q) 2 ): 

X'-p' l + l' 



H x > 
K" 
WJ 



r 



Hi + —K - I6ir— {finVn + X pV P ) , 
r r 



Ho , 1(1 + 1) 



r 

p A/2 



2r 



v> l + l 
~2 ~ 



e A/2 

X_ 8?r [imW n +xpW p ] , 

r 



, Ho + e x ' 2 rK-e x ' 2l ±±^V n -( l ^±+ n -)w n + 
I r \ r n J 



n 2 V% 



A (> 
npV% 

e x / 2 r 



e (A-«/)/2 rXn + n , (pv nWn + A° pW p ) 

e (X-u)/2 rXp + p , ( A nWn + C °pWp) 



H + e x/2 rK - e 



-v P - 



l + l P ' 
+ — 

r p 



W v 



P 2 Vl 
A 

npVl 



e {^-u)/2 rXp + p , ( C pWp + A o nWn 
(x -^' 2 rX n + n' (A° oP W p + B°nW n ) 



X n — 



I e y / 2 
r 2 



Ti/ji-i/) -n' (B^n + Alp) 



H a 



+nn 
+e»' 2 



r 2 
■' 1 



1(1 + 1) | ^ rc _ u/ 
' v 



Hi 



»n[---)-{B Q Q n + Alp) n' 



K 



ri (B Q nV n + A^pVp 



3 ( A -")/ 2 — n (BnW n + ApW p ) - Ane^+^Z 2 ^ (pnW n + XPW P ) 



-^(B°'nW n + A°'pW P )+ 



2n' A' - 



2r 



(23) 



(24) 



(25) 

(26) 
(27) 



(28) 



(29) 
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(B° nW n +A° oP W p )] 



(30) 



I 

— x p 

r 

+XP 
+e»' 2 



a u/2 



v/2 



Px(--v')-p'(C° Q P + A n) 



H 



1(1 + 1) 



- re 



-v/2 



Hi 



K 



fc^e^V (C° pV p + A° nV n ) 



Z 

z {x - u)/2 —p (CpW p + AnW n ) - 47re( A+ ^/ 2 ^ ( X pW p + [inW n ) 



,-{X-u)/2 



(C° P W p + A° nW n )] 



P_ 

r 



(cg'pWp + Jft'nWn) 



+ [ — + — p — 

V r ir r 



(31) 



The equations for the ordinary fluid envelope are obtained by taking the n — limit of the field equations. It is, 
however, not quite as straightforward as letting n — > in the above set of two-fluid perturbation equations. The 
reason for this is that we have in places divided through by which vanishes in the one- fluid limit. 

In the one-fluid case, the constraint equation becomes 



2-1 -I 2 



1 



o X-v 2 i a' 2 + ' ^ 

-2e u> + e - 



+lQTte x - l,/2 Xp = 



8tt* 



-,2A 



H - 
I - e~ A 



2uo 2 1(1 + 1) x fl- < 
~ ' 2 6 



8vr* 



8tt* 



I (l — e~ A ) - 1-kt 2 ^ 

2 v 1 



Hi 
K 



(32) 



The other two equations for the metric are 

a x 



H t ' 



-H 







A' — v 1 l + l 



k> = ^+ l M^ Hl 



Hi + 
l + l 



r 2r 

The final two equations are for the fluid and they take the form 

1 + 1 



e e 
— K — 167r — X pV, 
r r 

e A/2 

K - 8vr xpWj 



p ' 



Wp 



X p ' 



e X/2 r x/9 , /0 l(l + l) 

— — H + e x/2 rK - e x ' 2 y ' V p 



p ■ 



e (X-u)/2 r 

2^0 X P ' 



-—X„ 



v/2 



, 1 , 

PX\ v 

r 



pC p 



Ho +XP 



P u o 
e u ' 2 1(1 + 1) 



U -v/2 

— re 1 



Hi 



XP 



1 

2~i- 



n0 i 

C a pp 



K 



l -^e^ 2 P 'C° P Vp 



_ e (X-v)/2^ Cp 2 w _ 47re (A+,)/2 M!w 

r r 



+e -(X-v)/2 



-P-d'pWp+[^ 



2p> 



X' - v ' 
2r 



p' - ?-) C° oP Wp 



(33) 



(34) 



where 



X p =p 



W2 



X Ho + e- v l 2 uJ 2 (CpV p ) 



_ e {u-X)/2P_ (C*pW p ) 



(35) 



This set of equations is identical to that previously used |35| , p6| to study the oscillations of normal fluid neutron stars 
for a range of supranuclear equations of state. 

At the center of the star, the conditions are those given in Appendix A of Comer, Langlois, and Lin ^A, i.e. all 
functions are regular. At the surface, the conditions are the one-fluid conditions used by Detweiler and Lindblom 
p5| , |36| . The main difference here concerns the interface. The detailed treatment of the interface is discussed in 
Appendix I. It is found that the relativistic junction conditions imply that the three metric perturbations Ho, Hi and 
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K must be continuous at the interface. We assume that W n (R c ) is free to vary at the interface, i.e., the value it takes 
at the interface is determined by the general solution produced for the core. We also assume that X n (R c ) = 0, which 
will be shown below to be consistent with the chosen equation of state. From the results presented in Appendix I, 
it is found that two conditions remain. One of these implies that X p must be continuous at the interface. Recalling 
that we assume the proton number density and proton chemical potential to be continuous at the interface, i.e. that 
we have a second order phase-transition, the last condition is that W p must also be continuous at the interface. 



D. Computational strategy 

Even though our strategy for integrating the perturbation equations is similar to that used by Comer, Langlois, 
and Lin [Q there are some subtleties associated with the presence of the one-fluid envelope. Hence, it is worthwhile 
outlining the approach we have taken to the problem. In the core the perturbation equations can be written as the 
matrix equation 

dY 

-T = Q Y , (36) 
ar 

where 

Y = {H u K,W n ,W p ,X n ,X p } (37) 

is an abstract six-dimensional vector field. The 6x6 matrix Q depends on I, u>, the background fields, and the 
various coefficients A, A®, etc. As was shown by Comer, Langlois, and Lin one need only specify the set of values 
{K(0), W n (0), Wp(0)} at the center of the star. The remaining variables, Hi(0), X„(0), and X p (0), then follow from 
the r — > limit of the perturbation equations. All of the second derivatives, i?('(0), K"(0), etc, are likewise determined 
by {K(0),W n (0),W p (0)}. This means that, in order to provide information that enables us to construct the general 
solution, an integration starting from the centre must generate three linearly independent solutions Yi, Y 2 , and Y 3 . 
The corresponding general solution can thus be written 



Y(r)= »Y,(r) , (38) 



i=l 



where the Cj (i — 1, 2, 3) are constants to be determined. 

In the envelope, the problem is equivalent to the standard one for a single fluid and our strategy is identical to that 
of Lindblom and Detweiler [B5L B6| . We write the perturbation equations as 



dY 

dr 



Q Y , (39) 

wr 
where 

Y = {H U K,W P ,X P } (40) 

and the matrix Q can be deduced from Eqns. (^2|)-(|34|). At the surface of the star our solution must satisfy the single 
condition X p (R) = (cf. J3f], |3(|). This means that we must generate three linearly independent solutions in the 
envelope and the general solution can therefore be written 



Y(r)=j>Y,(r). (41) 



At the core/envelope interface we must enforce the additional condition that X n (R c ) = 0. We also know, from 
the analysis in Appendix I, that the variables Hi, K, W p and X p should all be continuous across the interface. This 
means that only the function value W n (R c ) remains unspecified. This is, however, as it should be since the interface 
represents a "free boundary" for the superfluid neutrons. (An analogy can be made with the water/air interface of the 
oceans, where the water oscillation is free at the interface.) Our strategy for solving the problem is based on two steps. 
First we continue the three solutions from the envelope into the core assuming in addition that W n (R c ) — 0. Then we 



9 



determine a fourth solution (needed to generate the general solution in the core) by assuming that W n (R c ) 7^ 0, but 
that all the other variables vanish at the interface. This means that we have determined a general solution of form 

7 

Y(r)=5>Yi(r) , (42) 

i=4 

where Y t (R c ) = Y t (R c ) for i = 4 - 6. 

The remaining step is to match the various solutions at some point r = r m in the core < r m < R c . At r m we 
must have 

3 7 

^c,Y 4 (r m ) = 5^CiYi(r m ) • (43) 

z=l i=4 

Once we have provided the overall normalisation (by specifying one of the c 4 coefficients), this problem can readily 
be solved for the remaining coefficients. This completes the solution of the interior problem. 

To determine the global solution we must also solve for the exterior metric perturbations. This problem reduces to 
that of integrating the Zerilli equation, cf. p2fl . Finally, in order to find a quasinormal mode (QNM) of the system we 
need to identify solutions that correspond to purely outgoing waves at spatial infinity. Our method for determining 
long-lived QNMs is described in Appendix III. 



III. AN ANALYTICAL ENTRAINMENT MODEL 

In order to extend the previous work by Comer, Langlois, and Lin further we want to construct a sensible 
equation of state that incorporates the entrainment effect. This effect arises whenever there is a coupling between 
two interpenetrating fluids, and has the net result that the momentum of one fluid is not simply proportional to that 
fluid's velocity. Rather, it is a linear combination of the velocities of both fluids. Hence, when one constituent starts 
to flow it will necessarily induce a momentum in the other. However, the entrainment effect is poorly understood, 
and there are as yet no completed relativistic models that we can base our discussion on. (The problem is currently 
being considered by Comer and Joynt, and we are hopeful that a fully relativistic formalism will soon be available.) 
The involved microphysics is, in fact, so uncertain that the best strategy corresponds to introducing some convenient 
parameterisation and then studying whether a variation in the chosen parameters affects the overall properties of the 
star and/or it's modes of pulsation. 

We noted earlier, in Sec. that the two chemical potential covectors /i^ and x M make manifest the entrainment 
effect by being linear combinations of the conserved four-currents, and p^ . Ultimately, this traces back to the master 
function A depending explicitly on the "entrainment variable" x 2 . Its presence in the background and perturbation 
equations listed in Sec. || is most apparent through the values of the coefficients A and A®. In the following we will 
use an expansion in terms of x 2 as outlined in Appendix II. At the heart of this expansion are the dimensionlcss 
ratios of the neutron and proton three- velocities with respect to the speed of light. This expansion makes sense 
since one would expect these ratios to be extremely small under most circumstances. Then the expansion becomes 
particularly accurate. The equation of state so constructed should have applications to the studies of quasinormal 
modes of neutron stars, with both static spherically symmetric and slowly rotating backgrounds. We thus expand the 
master function as 

00 

A(n^,x 2 ) = X>(rc 2 ,p 2 )(z 2 -^r ■ ( 44 ) 

where x 2 — np is expected to be small with respect to up. It should be noted that there are combinations other than 
x 2 — np that could be used in the expansion, the most obvious being those that have no dimensions, say, (x 2 — np) /np. 
These combinations are, however, not convenient in that the A, B, etc coefficients require that partial derivatives with 
respect to n 2 and p 2 be taken. The effect of this is that the expansions for A, B, etc would then not take the same 
form as Eq. ( fl4| ) above, since every derivative would bring in extra factors of x 2 outside the terms [(x 2 — np)/np\ % . A 
quick glance at Appendix II will show that if we use an expansion based on x 2 — np the basic form of the expansion 
is preserved for all coefficients that enter the field equations. 

In order to make contact with previous Newtonian studies of oscillating superfiuid stars we will adopt the particular 
entrainment model used by Lindblom and Mendell pi). To do this, a point of connection must be made between 
the coefficients used in the general relativistic superfiuid equations and the corresponding Newtonian equations. This 
connection can be made by taking the Newtonian limit of the general relativistic superfiuid formalism and then 
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comparing the dynamical variables (as in 0). In this way the mass matrix components p nn , p ppn and p np of the 
Newtonian formalism (cf. j37|) can be directly connected to the analytical entrainment coefficient Ai of the general 
relativistic formalism (cf. Appendix II) and the particle number densities n and p. Following this strategy we find 
that Ai can be written as 

c 2 m n m p 

Ai = 2~ —Pnp , (45) 

Pnp PnnPpp 

where 

m n n = p nn + p np , m p p = p pp + p np . (46) 

where m n r p \ is the neutron (proton) mass. The particular model of Lindblom and Mendell sets 

Pnp = -em n n , (47) 

where e is a constant. In the following we will refer to e as the "entrainment parameter." From the analysis of Prix, 
Comer and Andersson 1081 one can infer that 



m p p 



(48) 



where m* is the proton effective mass. Given that 0.3 < m*/m p < 0.8 (see, for instance, Sjoberg |3^]) "typical" values 
for a neutron star core may lie in the range 0.04 < e < 0.2. In the following we take this range as being "physically 
reasonable" . 

At this point it is relevant to note that, after all the numerical work described in this paper had been completed, 
Prix et al |f38| developed an alternative description of entrainment. While we could, in principle, have adopted our 
formulas here to this new description we have decided not to do this. Such a change would not have affected our 
discussion or the implications of our results in any way. 

IV. NUMERICAL RESULTS 

A. An analytical equation of state 

We will now apply our formalism to a simple model equation of state. We consider the case where each core fluid 
(at the lowest order in the expansion presented in Appendix II) behaves as a relativistic polytrope, i.e. we take 

A (rt 2 ,p 2 ) = -m„n - er, i n /3 " - m p p - apP^ . (49) 

This master function is clearly separable in n and p. Using the formalism developed in Appendix II, the relevant 
coefficients are (for the equilibrium configurations where x 2 = np) 

A = e 7" mp , ^ = 0, (50) 

m p p + e(m n n + m p p) 

B=^+a n p n n^ 2 -e m " m **/" , gg = aA {(3n _ i) n ^-2 (51) 

n m p p + e(m n n + m p p) 

C = ^ + a p P p p^ 2 -e ^m p n/p = _ /p _ 2 

p m p p + e(m n n + m p p) 

The pressure and chemical potentials of the core are 

* = <Jn (A, - 1) n - + a p (f3 p - 1) pft- , 

p = m n + (i n f3 n n p ' 1 ^ 1 , 



and 
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TABLE I: Parameters describing our stellar models I and II. Model I is identical to model 2 of j22j, and has no envelope. 
Model II, on the other hand, has an envelope of roughly 1 km and could be seen as a slightly more realistic neutron star model. 





model I 


model II 


u n /m n 


0.2 


0.22 


o~ P /m n 


0.5 


1.95 




2.5 


2.01 


Pv 


2.0 


2.38 


n c (fm -3 ) 


1.3 


1.21 


Pc (fm" 3 ) 


0.741 


0.22 


M/Mq 


1.355 


1.37 


R (km) 


7.92 


10.19 


Rc (km) 




8.90 



X = trip + cT p f3 p p " 1 . (53) 

Furthermore, one can verify that these coefficients and thermodynamic variables are such that the variable X n must 
vanish at the core/envelope interface, cf. Eqn. (p3|). 

Of course, we must also consider what happens in the envelope. This is much simplified because there cntrainment 
is not relevant. We get 

A(p 2 ) = -m p p - . (54) 

The relevant coefficients are 

C = ^+a pf 3 pP ^- 2 , C° Q =a p (3 p ((3 p -l) P ^ 2 , (55) 

and the pressure and chemical potential are given by 

* = a p {(i p -l)p^ , 

X = m p + (TpPpP "^ 1 . (56) 

We have considered two stellar models. The first model is identical to model 2 of Comer, Langlois, and Lin [f22| , 
and has no envelope. The second model has a significant envelope, and a more realistic (slightly smaller) proton 
fraction in the core. The relevant parameters that determine the two models are listed in Table |. In particular, the 
features of the core/envelope model are illustrated in Fig. [l]. The figure shows the radial profiles of the background 
neutron and proton particle number densities, n(r) and p(r), respectively. The model has been constructed to have 
the following features that are reckoned to be characteristic of real neutron stars: (i) A mass of about 1.4M Q , (ii) a 
total radius of about 10 km, (iii) an envelope of roughly one kilometer thickness, and (iv) a central proton fraction of 
about 10%. 



B. Quasinormal modes 

We have calculated the relevant fluid QNMs for our two stellar models, and various values of the entrainment 
parameter e. To put the results in context, it is useful to recall the basic features of the mode spectrum of a non- 
rotating, ordinary fluid neutron star. Despite the ordinary fluid system being the simplest possible description of 
the matter in a neutron star, there is nevertheless an impressive array of modes. In addition to the expected f- 
and p-modes, which are acoustic in nature, there also exists the so-called w-modes Q which are primarily due to 
oscillations of spacctime itself, with little coupling to the fluid of the star. If the equation of state has more than one 
parameter, and can be considered stratified — for instance, because of a varying proton fraction — then the star can 
also support low-frequency g-modes . 

One might expect that the essential doubling of the matter degrees of freedom due to the presence of a superfluid 
component might simply lead to a doubling of the families of modes of oscillation. However, it should not come as 
a great surprise that the w-modes in a superfluid neutron star look very much like those of the ordinary fluid case 
J22J. There is no doubling of modes: The w-modes are a feature of spacetime itself and depend on the curvature 
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FIG. 1: The radial profiles of the neutron and proton background particle number densities, n andp, respectively, for model II. 
The model has been constructed such that it accords well with a IAMq neutron star determined using the modern equation of 
state calculated by Akmal, Pandharipande and Ravenhall KG]. For reference, we show as horizontal lines the number densities 
at which Akmal et al suggest that i) neutron drip occurs, nj there is an equal number of nuclei and neutron gas, and iii) the 
crust/core interface is located. It should be noted that the latter should not coincide with our core/envelope interface since 
one would expect there to be a region where crust nuclei are penetrated by a neutron superfluid. 



induced by the background fluid rather than the actual nature of the fluid. But it is perhaps surprising that the simple 
expectation of mode doubling is not completely realized for the modes that are mainly due to matter oscillations. As 
mentioned in the Introduction, it has long been known that the additional fluid degree of freedom leads to the presence 
of a new set of modes, that have been dubbed superfluid modes. They are analogous to the ordinary fluid p-modes in 
that they are predominately acoustic in nature, (cf. Eq. (|l|)). The situation regarding g- modes is more confusing. Not 
only are there no new set of pulsating g-modes due to superfluidity, the standard set that one might expect to exist 
because of the varying proton fraction do not exist either. Inspired by Lee's pQ] numerical results, Andersson and 
Comer |DJ have used a local analysis of the mode spectrum to show that the g-modes disappear from the spectrum 
of pulsating modes. They prove that the "missing" g-modes exist as two independent sets of degenerate modes in the 
space of time-independent perturbations (the zero-frequency subspace). Their degeneracy will presumably be broken 
if one were to consider a two-fluid system governed by a three parameter equation of state. We plan to discuss this 
issue in detail elsewhere. 




-4 1 1— I 1 1 i 1 1 1 1 1 1 1 

0.1 0.2 0.3 0.4 0.5 

coM 

FIG. 2: This figure shows the asymptotic amplitude Aj n as a function of the (real) frequency uiM for our model I. The slowly- 
damped QNMs of the star show up as zeros of A[ n , i.e. deep minima in the figure. The first few "ordinary" and "superfluid" 
modes are identified in the figure, cf. Table III. 
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TABLE II: The oscillation frequency and the associated damping rate (in terms of the real and imaginary parts of uiM) for the 
first few fluid pulsation modes of model I. The w-modes for this model can be found in Table III of [22j. 



mode 


f 


so 


pi 


Sl 


J>2 


S2 


Re ljM 


0.137 


0.157 


0.306 


0.354 


0.585 


0.688 


Im toM 


7.2 x 10" 5 


4.2 x 10" 7 


6.6 x 10" 6 


5.2 x 10" 6 


3.2 x 10" 7 


5.5 x 10" 8 



Since the w-modes are largely unaffected by superfluidity, and we do not expect pulsating g- modes, we focus our 
discussion on the ordinary fluid f- and p-modes as well as the superfluid s-modes. We mentioned in the Introduction 
that there are two characteristics that distinguish the s-modes from their ordinary fluid counterparts: (i) counter- 
motion of the neutrons with respect to the protons, and (ii) a strong dependence on the parameters of entrainment. 
We will leave until the next subsection the discussion on the effects of entrainment. Hence, we first consider the 
spectrum of fluid modes for a model with A = 0. In Fig. || we provide a graph of the incoming wave amplitude 
(see Appendix III) versus the real part of the mode frequency for our superfluid model I. As discussed in detail in 
Appendix III, a QNM corresponds to those particular solutions where there is no incoming wave at infinity. Hence, 
the QNM frequencies correspond to the deep minima that can be seen in Fig. ^|. The main feature to notice is the 
presence of twice as many slowly damped fluid modes as in the ordinary fluid case, cf. Fig. 6 of p2| . 

In the earlier work of Comer, Langlois, and Lin there was no attempt at determining the gravitational- 

wave damping times for either the ordinary fluid or the superfluid modes. The reason for this was the difficulty in 
determining complex QNM frequencies with imaginary parts that are orders of magnitude smaller than the real parts. 
To deal with this problem, we have developed a new method, which is outlined in Appendix III. The results we have 
obtained for model I are given in Table O. As one can see, the ordinary fluid modes have damping rates that are 
consistent with what is known from calculations of ordinary fluid neutron stars |3j], ||(| . It should also be noticed 
that most of the superfluid modes have gravitational-wave damping rates that are similar to the high order p-modes. 
This is to be contrasted with recent statements in Ref. 



fact, in Section [v| below we will provide a proof that all ( 



43 that superfluid modes will not radiate gravitationally. In 



1 QNMs of a superfluid star must radiate. This means that 
the superfluid modes could, at least in principle, be relevant for gravitational-wave asteroseismology. This possibility 
will be discussed further in Section VI. 

So far we have mainly discussed the results obtained for model I, for which there is no envelope. When we turn to 
model II, which has an ordinary fluid envelope of roughly one km, we find that the results do not change qualitatively. 
In particular, we do not find that new modes arise because of the presence of the envelope. At first this may seem a 
little bit surprising. Especially since it is well known that an elastic crust supports several additional sets of modes. 
However, in our case the absence of new modes is due to the nature of the phase transition at the core/envelope 
interface. We have chosen the phase-transition to be second order, which means that the number density of the 
protons remains smooth as the superfluid neutron component vanishes (at r — > R c ). Should we have taken the phase- 
transition to be first order, i.e. allowed for a jump in the proton number density at R c , our calculations would have 
unveiled a set of interface g-modes. 

C. The effect of entrainment — avoided mode crossings 

The two main goals of the work presented in this paper were: i) to allow for the presence of a core/envelope 
transition, and thus in principle be able to consider cases where the superfluid constituent is confined only to a part 
of the star, and ii) to determine how entrainment affects the QNM frequencies. As we will now discuss, the effect of 
entrainment can be considerable. 

We have carried out a series of calculations for model II and entrainment parameters in the "physical range" 



0.04 < e < 0.2. A sample of results are given in Table III. Listed there are the oscillation frequencies and associated 
damping rates for the first few pulsation modes of model II and three different values of e. First of all it is relevant 
to compare the results for the case of vanishing entrainment to those for model I, cf. Table [n|. Recall that the main 
difference between models I and II is that the latter includes an ordinary fluid envelope. Nevertheless, it is clear from 
the numerical results that the QNMs of the two models are qualitatively quite similar. When we turn to the effect of 
varying the entrainment parameter we find that the superfluid mode frequencies shift considerable, while the ordinary 
fluid modes remain virtually unchanged. This is not a surprising result given the discussion in Ref. ]i"7| and Eq. (Q). 
It is also relevant to note that the gravitational-wave damping rates can be strongly affected by a change in e. 

As we will now discuss, the effect that a varying entrainment has on the damping rate of the modes can be 
understood from the results illustrated in Figs. || and |J. Fig. || shows how the oscillation frequencies for the first few 
modes change as the entrainment parameter is varied within the acceptable range. The modes that have ordinary fluid 
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TABLE III: The oscillation frequency and the associated damping rate (in terms of the real and imaginary parts of ujM) for 
the first few fluid pulsation modes of model II. We show results for three different values of the entrainment parameter e. These 
correspond to the case of no entrainment as well as the upper and lower limits for the range that we take as "physically realistic" . 
From this data we see that, while the ordinary fluid modes are hardly at all affected by the entrainment, the superfiuid mode 
frequencies vary by as much as 10% as e is varied within the realistic range. 





mode 


f 


so 


Pi 


Si 


P2 


S2 


P3 


S3 


€ = 


Re ujM 


0.112 


0.124 


0.252 


0.288 


0.379 


0.424 


0.502 


0.554 




Im ujM 


4.9 x 10 -5 


8.9 x 10 -6 


1.3 x 10 -5 


1.4 x 10 -5 


8.0 x 10 -8 


1.3 x 10 -7 


1.5 x 10 -8 


1.7 x 10 -9 


e = 0.04 


Re ujM 


0.113 


0.130 


0.253 


0.299 


0.382 


0.442 


0.507 


0.577 




Im ujM 


5.3 x 10" 5 


5.3 x 10" 6 


1.4 x 10" 5 


8.5 x 10" 7 


9.2 x 10" 8 


1.2 x 10" 7 


1.2 x 10" 8 


4.3 x 10" 9 


e = 0.2 


Re cuM 


0.113 


0.149 


0.257 


0.328 


0.394 


0.491 


0.523 


0.639 




Im ujM 


5.5 x 10 -5 


3.2 x 10 -6 


1.4 x 10 -5 


3.2 x 10 -7 


1.0 x 10 -7 


1.2 x 10 -7 


2.3 x 10 -9 


1.2 x 10 -8 



behavior, i.e. for which the neutrons and protons "flow together" , in the limit of vanishing entrainment arc shown as 
solid lines in the figure whereas the superfiuid modes, where the neutrons and protons are largely counter-moving, are 
given by the dashed lines. The main feature of the figure is the presence of so-called avoided crossings. For the higher- 
order modes (near the top of the figure) there are points in the (e, Re ujM) plane where the solid and dashed lines 
approach each other, but rather than crossing they diverge from each other. An interesting aspect of these avoided 
crossings can be gleaned from Fig. ^, which shows the Lagrangian variations An and Ap for the particular modes that 
correspond to the points a e and b e , e = 0, 0.1, 0.2, of Fig. |^. For the mode a the two fluids are largely counter-moving 
as e — > 0, and for mode b the two fluids are then comoving. However, as e -> 0.2 we see that the modes have changed 
character in that it is now mode a that is comoving and mode b that is counter-moving. This exchange of character 
of modes is a characteristic of avoided crossings familiar from other problems in stellar pulsation theory E3|. The 
presence of avoided crossings between stellar pulsation modes is familiar from many other situations, but we believe 
that our results provide the first hard evidence for the presence of this phenomenon in superfiuid stars. 

From general post-Newtonian arguments, one would expect the co-moving modes to radiate gravitational waves 
more efficiently than the counter-moving superfiuid modes. Thus it is not surprising to find that the damping rate of a 
superfiuid mode increases as entrainment is varied and an avoided crossing is approached. In the present context, this 
effect is probably not distinct enough to be of great relevance but one can argue that it may be of great importance in 
closely related problems. As has been argued recently [|17j , avoided crossings may be at the heart of recent calculations 
of the effect of superfluidity on the r-mode instability. Lindblom and Mendell ]f4j have analyzed the effect of mutual 
friction damping on the r-modes using the same entrainment model as we employ here. They found that mutual 
friction was, in general, not effective at damping the r-modes. However, they also found (cf. their Fig. 6) that the 
mutual friction damping time could be very small for particular values of e. We believe that a proper explanation of 
this result can be obtained via the avoided crossings phenomenon. The basic idea is that mutual friction should be 
most effective whenever the neutrons and protons are counter-moving, as with the superfiuid modes. Andersson and 
Comer |L7| have shown that there are two sets of r-modes, that are quite analogous to the ordinary fluid and superfiuid 
modes discussed in this paper in the sense that one set has the neutrons and protons comoving whereas the other set 
has them countermoving. Although avoided crossings between these two classes of r-modes in a superfiuid star have 
not yet been demonstrated, we believe that our current results provide strong support for the idea by demonstrating 
the presence of avoided crossings in a very closely connected situation. We thus assert that the particular values of 
e for with Lindblom and Mendell find small mutual friction damping times correspond to stellar models for which 
the two classes of superfiuid r-modes are close to an avoided crossing. The veracity of this argument remains to be 
confirmed by detailed calculations that we plan to carry out in the near future. 

V. ARE THERE NON-RADIATIVE MODES IN SUPERFLUID NEUTRON STARS? 

In this Section we digress somewhat and focus our attention on a question of principle: Is it possible to have QNMs 
that do not radiate gravitationally in a superfiuid star? The question is motivated by the simple fact that, while every 
non-radial motion induced in a one-fluid system must lead to the emission of gravitational waves, the situation could 
conceivably be different in the two-fluid case. One can imagine the possibility that the two fluids move in such a way 
that the average mass-density flux vanishes identically. In fact, as we have discussed elsewhere jl^], the superfiuid 
modes are essentially of this nature and it has been argued by Sedrakian and Wasserman 0] that this means that 
these modes will not radiate. Of course, one must not forget that, in the post-Newtonian picture, gravitational waves 
are associated with both mass- and current multipoles. Thus even the extreme case of a two-fluid oscillation such 
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FIG. 3: This figure shows how the frequencies of the fluid pulsation modes for our model II vary with the entrainment 
parameter e. The modes shown as solid lines are such that the two fluids are essentially comoving in the e — > limit, while 
the modes shown as dashed lines are countermoving. As is apparent from the data, the higher order modes exhibit avoided 
crossings as e varies. Recall that the range often taken as "physically relevant" is 0.04 < e < 0.2. We indicate by a e and b t the 
particular modes for which the eigenfunctions are shown in Fig. W. 



that the mass-multipoles vanish identically is likely to radiate through the induced currents. Thus it may be very 
difficult to set the two fluids into motion without generating gravitational waves. And this is, indeed, the way that it 
turns out. As we will prove below, a two-fluid star can (in general) have no non-radiative modes. 

We consider the superfluid perturbation equation in the special case where the metric is left completely unperturbed, 
i.e. when no gravitational radiation is created. Assuming a given, fixed non-rotating background, we can obtain the 
corresponding field equations by setting to zero the three metric perturbation Hq, Hi, and K. This results in the 
following nine equations for six matter variables, cf. Eqns (^5[)-(|3l|), 
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FIG. 4: An illustration of the fact that the modes exchange properties during an avoided crossing. We consider two modes, 
labelled by a e and b t (cf. Figure []) . The mode eigenfunctions are represented by the two Lagrangian number density variations, 
An and Ap (solid and dashed lines, respectively). For mode a the two fluids are essentially countermoving in the e —* limit (it 
is a superfluid mode), while the two fluids comove for mode b (it is similar to a standard p-mode). After the avoided crossing 
(which takes place roughly at e — 0.1) the two modes have exchanged properties. 
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(C°pW p + A° nW n 



X n = ne-^ 2 Lo 2 (BnV n + ApV p ) - c>- A )/ 2 - (B° Q nW n + A° pW p ) 



X p = pe-" /2 Lu 2 (CpVp + AnV n ) - e (l/ - A)/2 - (C°pW p + A° nW n ) 



(62) 



(63) 



(64) 



(65) 



It is clear that, unless some of these equations are linearly dependent the problem is overdetermined and we cannot 
have a non-trivial solution. In order to show that the problem is in general overdetermined, we use the first three 
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equations above to rewrite W pi V p , and X p in terms of the neutron variables W n , V n , and X n , and then substitute 
these into the remaining six equations. Remarkably, one finds that Eqs. ( |62[ ) and ( |63| ) yield exactly the same result 
when this substitution is done. One also finds that Eqs. (^J) an d ( Pj ) yield the same result after the substitution. 
We are thus left with Eqs. ( |60| ) and (|6l|). After the substitution, these equations yield the following two equations for 
W n : 



_ e v,Mi±i) Vn _ i±V„ + <U. Ed] - 1 

r r VI V X J \n p 



r T>" V 

— 2 7 ) 

In the manipulations we have used repeatedly the background equations (|ll|) and 

M' = "\^ v ' > x' = -^X^' • (68) 

In order to complete the argument in a clear way, we restrict ourselves to the case of vanishing entrainment and 
separable equations of state, i.e. we concentrate on the case A = A% = 0. This is not at all necessary, but it simplifies 
the analysis considerably since the involved relations are much simpler. Taking the difference between ([36|) and d67j ) 
we require 

Bgcg (- - P - V„ + (c°p - *S °n) X„ = (69) 

\ n p ) pn 2 \ P J 

in order for the problem not to be over-determined. Using the definition of X n we can rewrite this equation as 

=V2-^,,,2 R ( r _ !!X R o\ v _ , ^0 („t„0X _ Wr 



e^ru?B C ° - ^B° )V„ + — [ n'B° ± - p'C° Q )W„ = 0. (70) 
V PM / P \ P J 

Now using the fact that x = and /i = Z?ri in the case of vanishing entrainment, we have 

e A / 2 -^ 2 (£C ° - CB° ) V n + ^ (C£ ^ - 6C ^ TU„ = e^roj 2 {BC° Q CB° ) V n = . (71) 



To arrive at the last equality we have employed the background identities (|ll|) again. Clearly, we are now left with 
two possibilities. Either V n vanishes, or the equation of state must be such that 

BC° - CB% = . (72) 

In the first case one can show that the corresponding solution for W n is 

W n <x — [xt . (73) 
nr ~*~ 

This solution is clearly physically unacceptable since it diverges at the centre of the star. In other words, if V n = 
we must also have W n = 0, i.e. the trivial solution. It is easy to show that the second case requires that the master 
function be such that 

g a dk d / dA \ _ 8n2 dA__d_ f dh_ \ _ Q 
dn 2 dp 2 \ dp 2 J dp 2 dn 2 \ dn 2 J 

This is clearly a very particular form. We have thus proven that the QNMs of a superfluid star must radiate unless 
the equation of state belongs to a very special class. The obvious exceptions are i) when the master function takes 
the form 

A = <r n n 2 + (j p p 2 (75) 
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and ii) whenever A depends on n and p in identical ways. 

It is worth pointing out that the calculation in this section was carried out within the Regge- Wheeler gauge, and 
given possible gauge issues one may worry that this means that the result is of limited validity. However, the result 
will hold in general since we can easily construct gauge- invariant quantities (such as the Zerilli function) that represent 
the gravitational-wave degrees of freedom from the metric perturbations calculated in any particular gauge. In our 
case it is trivial to see that, if all metric perturbations vanish identically the Zerilli function will be identically zero 
and no gravitational-waves will emerge from the system. 



VI. DETECTABLE GRAVITATIONAL WAVE SIGNALS? 



Given that a new generation of gravitational-wave detectors are likely to be operating at their projected levels 
sensitivity within the next few months, it is appropriate to conclude this paper with a brief discussion of a possible 
future application of our results. Suppose that the various modes of a superfiuid neutron star were excited to an 
amplitude such that the associated gravitational-wave signal could be detected. To what extent would it then be 
possible to solve the inverse problem and deduce information concerning the superfiuid parameters? In other words, 
can we hope to use "gravitational-wave asteroseismology" to probe the superfiuid interior of mature neutron stars? 
This question has recently been discussed by two of us p^ j , and therefore we only provide a brief background here. 

As has already been discussed by Kokkotas et al S , the detection of gravitational- wave signals from pulsations in 
newly born neutron stars could be used to infer the mass and radius of the star. This information would put strong 
constraints on the supranuclear equation of state. However, this argument relies on releasing an energy equivalent 
to something like 10~ 5 Mqc 2 through the QNMs. This may not be unreasonable for the wildly pulsating object 
formed through a strongly asymmetric supernova collapse, but it is difficult to think of a mechanism whereby the 
oscillations of a mature (and thus superfiuid) neutron star core will be excited to a similar level. Instead, we take as 
a "reasonable" order of magnitude estimate the energy associated with a typical pulsar glitch. The released energy 
can then be estimated as 



AE re IQAQ. « (1CT 6 - 10~ 8 )JO 2 



(76) 



where ft — 2ir/P is the rotation rate of the star, and P is the observed pulsar period. In this formula it is appropriate 
to use the moment of inertia / ~ 10 45 gem 2 of the entire star, since the spin-up incurred during the glitch remains 
on timescales that are much longer than the estimated coupling timescale between the crust and the core fluid. By 
combining the above formula with the data for typical glitches in the Crab and Vela pulsars, cf. Table IV, we arrive 
at estimates of the energy associated with typical glitches that accord well with suggestions in the literature Gsl pi . 



TABLE IV: Data for archetypal glitching pulsars. 



PSR 


P (ms) d (kpc) Afl/n 


AE/Mqc 2 


Crab 
Vela 


33 2 10" 8 
89 0.5 10" 6 


2 x 10" 13 

3 x 10~ 12 



As is evident from Table IV , we expect a glitch to be associated with energies of the order of 10 -13 — 10~ 12 M Q c 2 . In 
the following we will assume that a similar energy is channeled through the various pulsation modes. This then allows 
us to use the formulas obtained by Kokkotas et al j5j to estimate, for a given detector configuration, the attainable 
signal-to-noise ratio. 

Assume that the gravitational- wave signal from a neutron star pulsation mode takes the form of a damped sinusoidal, 
i.e. 



h(t) = Ae- {t ~ T)/td sin[2vr/(t - T)] for t > T 



(77) 



where / is the frequency of the QNM, td is its characteristic damping time, and T is the arrival time of the signal at 
the detector (thus h(t) =0 for t < T). Using standard results for the gravitational- wave flux H, the amplitude A of 
the signal can be expressed in terms of the total energy radiated through the mode: 



A & 7.6 x 10' 



AEq 1 s / 1 kpc 



1 kHz\ 



io- 12 U 

where AE Q = AE/M Q c 2 . The signal-to -noise ratio for this signal can be estimated from 

4Q 2 AH d 



1 + 4Q 2 2S n 



(78) 



(79) 
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TABLE V: The frequency and damping rate for the first few modes of our Model I (which is identical to model 2 of Comer 
et al We also show the gravitational- wave signal-to- noise ratios resulting from the "glitch model" discussed in the main 

text. The results correspond to an advanced EURO detector with (case 1) and without (case 2) photon shotnoise, respectively. 
The lower estimate is for a Crab glitch while the upper estimate follows from the Vela data. 



Mode 


f (kHz) U (s) 


Case 1 Case 2 


/ 
Pi 


3.29 0.092 
7.34 1.01 


0.4—6 300—4.7 x 10 3 
0.08—1.2 680—1 x 10 4 


so 

Sl 


3.76 15.75 
8.49 1.29 


0.3—4.8 350—5.4 x 10 3 
0.06—0.9 780—1.2 x 10 4 



where the "quality factor" is Q = nftd and S n is the spectral noise density of the detector. 

We can now combine these estimates with the QNMs from (say) Table |TJ. The relevant frequencies and damping 
times in kHz and ms, respectively, are given in Table 0. When we compare the obtained estimates to the new 
generation of large-scale interferometric detectors it immediately becomes clear that these signals would be too weak 
to be detected. This is illustrated in Figure || where we compare the dimensionless strain \ffS n for various detector 
configurations to the estimated strain caused by the QNM signals. 

This does not, however, mean that we should give up on the main idea behind this analysis. We simply have to 
acknowledge that we are likely to require significant improvements in technology if we are to be able to make this 
into a viable approach. But it seems inevitable that the available technology will improve over the next decades. In 
fact, various groups are already discussing possible improvements in detector sensitivity that may be achievable in the 
future. In order to illustrate the levels that are being discussed, we consider the so-called EURO detector, for which 
the sensitivity has been estimated by Sathyaprakash and Schutz (for further details see [|7) ) . We consider two possible 
configurations: In the first, the sensitivity at high frequencies is limited by the photon shot-noise, while the second 
configuration reaches beyond this limit by running several narrowbanded (cryogenic) interferometers as a "xylophone" . 
The corresponding noise-curves are illustrated, and compared to the current generation of interferometers, in Figure]^. 
It is immediately clear from Figure |^ that a EURO detector would be a superb instrument for studying pulsating 
neutron stars. This means that previously suggested strategies [|] for unveiling the supranuclear equation of state may 
eventually be put to the test. In fact, as is clear from the estimates in Table M, where we list the signal-to-noise ratio 
estimated from ((79|), one should also be able to detect the superfiuid oscillation modes. The various modes would be 
marginally detectable given this level of excitation and a third generation detector limited by the photon shotnoise. 
If this limit can be surpassed by configuring several narrowbanded interferometers as a xylophone, the achievable 
signal-to-noise ratio will be excellent. Hence, it seems plausible that one could infer the parameters of neutron star 
superfluidity. As an obvious "by product" one might hope to shed light on the mechanism for pulsar glitches. 




FIG. 5: The spectral noise density for the new generation of laser-interferometric gravitational-wave detectors that will come 
online in the next few years (thin lines) is compared to speculative estimates for the futuristic EURO detector (solid lines). A 
key feature of this advanced configuration is that it may operate several narrowbanded interferometers as a xylophone, thus 
reaching high sensitivity at kHz frequencies. We also indicate the effective gravitational- wave amplitudes from the glitch-induced 
mode-oscillations discussed in the text. 



We can also confirm that, provided that the modes can be detected, the oscillation frequencies can be extracted 
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with good accuracy from the data. By combining Eqns. (11) and (12) from Kokkotas et al || with ( |79| ) one can 
obtain a formula for the signal-to-noise required to determine the mode frequency with a relative error crj/f. We get 
(since Q >> 1 for the modes under consideration) a relation 

From this relation we can deduce that a detection with signal-to-noise of (say) 10 would enable one to infer the fluid 
mode frequencies with an accuracy of the order of a percent or so. With this level of precision one should be able to 
distinguish clearly between the "normal fluid" f and p-modes and the superfluid s-modes. In other words, we would 
have the information required not only to infer the mass and radius of the star Q , we could also hope to constrain 
the parameters of neutron star superfluidity. 

VII. CONCLUSIONS 

The main motivation behind the work presented in this paper is the fact that neutron star physics is not adequately 
modelled within Newtonian gravity. It is well known that Newtonian results differ greatly from the correct relativistic 
ones already at the level of determining the mass and radius of a star with a prescribed central density from a given 
"realistic" equation of state. As far as the QNMs of the star are concerned, Newtonian studies are extremely useful 
since they help us understand the physics of different classes of modes. But at the same time it is clear that if we 
require a detailed model of the oscillation spectrum of an astrophysical neutron star we must approach the problem 
from the relativistic point of view. This is particularly crucial if we are interested in the gravitational-wave damping 
rates. Furthermore, since mature neutron stars are likely to contain several superfluid components, it is important 
that we develop a framework for modelling multi-fluid systems in full general relativity. The present paper represents 
significant progress towards this goal. 

We have considered a core-envelope model, with superfluid neutrons being present only in the core of the star 
while the outer region is composed of an ordinary fluid. Our model includes a simple, yet reasonable, model for 
entrainment and we have shown how neighbouring QNMs undergo avoided crossings as the entrainment parameter is 
varied. Finally, we have ruled out the possibility that a two-fluid star could have non-radiative pulsation modes, and 
discussed the possible detection of gravitational-wave signals from oscillating superfluid neutron stars. 

In the near future, our aim is to turn our attention to the oscillations of rotating superfluid stars. This is an exciting 
problem area, since various modes of oscillation may be driven unstable by the emission of gravitational waves. Of 
particular recent interest are the so-called r-modes, and the damping due to superfluid mutual friction. In a recent 
work, Lindblom and Mendell [T4j showed that mutual friction was effective at damping unstable r-modes only for very 
particular values of the entrainment parameter. Andersson and Comer Jl7| have speculated that this peculiar behavior 
might be due to the existence of avoided crossings between two classes of r-modes (analogous to the two classes of 
acoustic fluid modes discussed in the present paper). The results of the present paper lend support to this idea. Yet 
it is clear that our understanding of the various involved issues is far from complete, and that a considerable amount 
of work remains to be done before we can model the dynamics of rotating superfluid neutron stars in a "realistic" 
way. 
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Appendix I: The Junction Conditions 

In this Appendix will be presented a geometrical approach to deriving the junction conditions that must be used 
to smoothly join together an inner superfluid core with an exterior normal fluid envelope. The junction conditions 
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will be obtained via an analysis of the first and second fundamental forms associated with the (in general, timelike) 
hypersurfaces given by the level surfaces of the generalized pressure "3/ (cf. Eq. (||)) If there are no "delta- function- 
like" discontinuities in the pressure, then the first and second fundamental forms will be continuous throughout the 
star [E5). Turning the problem around, by demanding the continuity of the first and second fundamental forms a 
smooth joining at the core/envelope interface can be achieved. 

Consider that the level surfaces of "J" are timelike, i.e. the normal to these surfaces is given by 

W 7 ' = 9 V " , (81) 

where Af^JV^ = 1. The so-called first fundamental form 7^ (i.e. the induced three-metric) of these level surfaces is 

7 „„ =±%±l g aT , (82) 

where the "perp" operator is given by 

^=^-^<W„. (83) 

The second fundamental form (i.e. the extrinsic curvature) of the level surfaces is defined as 

= - ±l±l V (a N T) , (84) 

where the parentheses imply symmetrization of the indices. 
Let us consider the pressure to be of form 

*(t,r,fl) = *o(r) +6*(t,r,0) . (85) 

The components of the unit normal are thus found to be 

W° = - e Va-^ + e -(A/a+«)^ 01 , 
r 2 % ' 

TV 3 = . (86) 
The non-zero components of the first fundamental form are 

700 = -e + dgoo , 701 = og i - e ^7- , 712 = -e , 

722 = r 2 + £.922, 733 = sin 2 6Y 2 + <5g 33 , (87) 

and we also find 

^ y (V /2 |j - , ^02 = [ ^£ - ^ 50M 

^ 12 = — ~W K22 = _ iv* e % " ^ ^ 922 " 

if 33 = -sin 2 6> -j^ + cot6»e A/2 — f + ^—Sc/ 33 - ——6g n 

y e A/2 ^ 2sin 2 6»e A / 2 2e 3A / 2 

for the non-trivial components of the second fundamental form. 

Given that a smooth background can be constructed independently of the oscillations, we can assume that the 
background and linearized pieces of 7 Mi , and K^ v are separately continuous at the core-envelope interface. We will 



-v- 








5^,9 


1 


% 
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consider the background pieces first. For clarity of presentation the matter and metric variables in the envelope will 
be distinguished by a tilde. The continuity of the first fundamental form thus yields for the background 



A(i? c 



v{R c ) = u(Rc) 
and continuity of the second fundamental form implies 

u'(R c ) = v'{R c ) 

Using the background Einstein equations, and defining 

e~ A = 1 - 2m(r)/r 

where 



e- x = l 



oKRc) 



2rh(r)/r 



m(r) 



-47T 



r 2 Ao(r)dr , m(r) 



-47T 



r 2 A {r)dr + C 



(89) 



(90) 



(91) 



(92) 



we find that the background junction conditions imply C = m(R c ) and ^q(R c ) = ^>o(R c ), where ^>q and — Ao are the 
pressure and energy density of the envelope, respectively. It is also useful to note that the Tolman-Oppcnhcimer-Volkov 
equations for the core and envelope imply 



%(R C ) V (R C )-A (R C 



%(R C ) y (R c ) - A {R C ) 



(93) 



Before dealing with the linear perturbations, it will be convenient to write out the linearized pressure as a function 
of the fundamental matter and metric variables, for both the core and the envelope. We have used the field equations 
to help simplify the formulas. The final forms are (for the radial dependence) 



5^ 



1 



^aTJ 0™W n + XPW P ) ~— 2 {X n + X p ) 



and 



5^> 



2re A/: 



;XPW P 



oV/2 



(94) 



(95) 



where we have again put a tilde over all the linearized variables associated with the envelope. Now, it can be seen 
that the junction conditions imply that the metric perturbations are continuous at the core/envelope interface, i.e. 

Hq(R c ) = Hq(R c ) , 
H^R.) = HiiRc) , 

K(R C ) = K(R C ) . (96) 
The matter variables, on the other hand, must satisfy two conditions, which are 

x(Rc)P{Rc)W p (R c ) = fi(R c )n{R c )W n (R c ) + X {Rc)p(R c )W p (R c ) (97) 

and 



X P (R C ) 



y {R c ) - A (R C ) 
*o(Ri) - A (R C ) 

x(Rc)P(Rc)W p (R 



(X n (R c ) + X P (R C )) 



vl{Rc)e (v{R c )-HR c ))/2 / y o{Rc) _ ~ Aq{Rc) 

2R C { MRo)- A Q (R C ) 



(98) 



It is important to notice that the junction conditions do not imply that 5^ = 5^> at the core/envelope interface, 
but rather that S^S/^q = S^f/^!' . This fact is the clearest demonstration why a geometrical approach is crucial 
for obtaining the correct junction conditions, since this particular result is a direct consequence of the fact that 
Schwarzschild-like coordinates cause derivative discontinuities whenever the energy density is not continuous. If it 
is the case that Ao(i? c ) = Ao(i? c ), then we can see from Eq. ( p3| ) that the junction conditions will imply 5^ = S^>. 
Having a continuous energy density also implies that X p (R c ) — X n (R c ) + X p (R c ). 



23 



Appendix II: The Analytical Equation of State 

In this Appendix we develop the analytical equation of state used in the main part of the text. The essential 
strategy is to introduce an expansion based on the assumption that the fluid velocities are small compared to the 
speed of light. This is a reasonable assumption for neutron star pulsations. 

A. A Local Analysis of the Entrainment Parameter 

Recall that the entrainment variable x 2 is given by x 2 = —n p p p . It is convenient to write each conserved four-current 
as in the main text: 

n p = nu p , p p = pv p , (99) 

except that now we take u p u p = —c 2 and v p v p = —c 2 where c is the speed of light. If r„ and t p denote the proper 
times of the neutron and proton fluid elements, respectively, then the worldlincs of each fluid element are obtained 
from the functions 

<(r„) = (i(r„),<(r„)) , x£(t p ) = (t(r p ) , x p (r p )) . (100) 
The "unit" four-velocities are thus given by 

dr^ dr^ 
u p = ^ , = ^ . (101) 

d,T n dTp 

Consider, for the moment, a region within the fluid that is small enough that the gravitational field does not change 
appreciably across the region. In this case, a locally Minkowski coordinate system can be set up and the metric can 
be approximated by the flat metric: 

ds 2 = -d{ct) 2 + 5 l0 dx l dx3 . (102) 

Letting u l n = dx l n /dt and v p = dx p /dt, as well as u 2 = SijU l n u J n and v 2 — 8ijV p v p , one can show that the four-velocity 
components can be written as 

u° = , C , u l = - , (103) 

y/l - K/c)2 v/1 - K/c) 2 

and 

v° = - C , v l = - V l . (104) 

With this decomposition one finally obtains for the entrainment variable 

If it is the case that the individual three-velocities are small with respect to the speed of light, i.e. that 

M«i , (106) 

c c 

then it will be true that x 2 sa np to leading order in the ratios u n /c and v p /c. This basic fact will be at the heart of 
the expansion considered below. 

B. The Analytical Equation of State 

Given what was just discussed, it makes sense to consider equations of state that can be expanded like 

oo 

A(n 2 ,p 2 ,x 2 )=J2^(ri 2 , P 2 )(x 2 -npy , (107) 

i=0 
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since x 2 — np is small with respect to np. In terms of this expansion, one can show the following for the A, Aq, 
etc. coefficients that appear in the field equations: 

oo 

A = \i(n 2 ,p 2 ) (x 2 ~ up) 1 

z=l 

B = -i^l-lA-^^i^-npy , 
n on n n — ' on 

i—i 

p op p p ^—^ op 

■ d 2 \ Q y—^ 9 2 \i , 2 si 

opon upon v ' 

2 — 1 

It is important to note that „4q vanishes if the master function is such that A^ are separable in n and p. 

The utility of this expansion is especially apparent for the quasinormal mode calculations, because when any of the 
coefficients are evaluated on the background, then one sets x 2 = np, and thus only the first few A^ are needed. In 
fact, in our analysis we only retain Ao and Ai, where the latter contains the information concerning the entrainment 
effect. 



Appendix III: An accurate method for determining long-lived QNMs 

The problem of calculating quasinormal modes of relativistic systems, such as black holes and neutron stars, is in 
many ways far from trivial. In general, a strategy for finding QNMs has to involve a prescription for imposing a pure 
outgoing wave boundary condition at infinity for a linear second-order differential equation. In the context of the 
present paper the relevant equation is the one derived by Zerilli : 



Z = . (109) 



Here r* is the standard tortoise coordinate and we have assumed that the perturbed quantities have a harmonic 
dependence on time, i.e. behave as exp(iu)t). The effective potential V(r) is rather complicated but here we need only 
know that it is such that the behaviour of a general solution at spatial infinity (as r* — ► +oo) is 

Z ~ A oat {u)e- iur ' + A hl (cu)e^ r ' . (110) 

A QNM of the system is a solution that combines some physical constraints (no waves coming out of the event 
horizon in the case of a black hole or a regular solution to the equations for the interior of a star) with purely outgoing 
waves, A- U1 (uj n ) — 0, at infinity. Two typical difficulties arise in the determination of such mode solutions. Both are 
due to the fact that the QNMs are damped by gravitational radiation emission, and thus the QNM frequency cu n 
must be comple x (wi th a positive imaginary part unless the mode is unstable). This means that an outgoing- wave 



solution to Eq. (109) will be exponentially growing as r* — > +oo and one would, in principle, need exponentially 
high numerical precision to discard the ingoing solution. This problem is particularly challenging for rapidly damped 
modes, like those of black holes and the neutron star w-modes. A second difficulty that arises is relevant for the 
p-modes of a neutron stars. The high order p-modes are damped very slowly by gravitational radiation. Thus, the 
characteristic frequencies are such that the imaginary part is many orders of magnitude smaller than the real part. 
Although not as conceptually challenging as the first QNM problem, the difficulty associated with extremely small 
imaginary parts typically prohibits the determination of any but the first few of the neutron star p-modes. In this 
Appendix, we describe a new method for dealing with this problem. The various mode results presented in the main 
body of the paper were obtained within this scheme. 

Let us assume that our relativistic system has a QNM with complex frequency oj n . Then the solution to the 
Zerilli equation is such that Ai n (tu n ) = as r* — > oo whereas Ao U t(u> n ) ^ asymptotically. From the fact that it is 
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the frequency squared that appears in all the relevant perturbation equations (see Sec. 0) we can draw two general 
conclusions: i) another outgoing- wave mode is characterized by —Cu n (where the bar represents complex conjugation), 
and ii) time-reversal of an outgoing- wave mode leads to a solution that corresponds to purely ingoing waves at infinity. 
That is, a solution such that, at infinity, A out (— w n ) — while Ai n (—uj n ) ^ 0. Moreover, it follows that a second such 
solution corresponds to u> n . As we will demonstrate below, this information can be very useful when trying to identify 
normal-mode frequencies that have small imaginary parts. 

Most schemes for determining QNMs are based on numerical constructions of the ratio k — A out /Ai n . Assuming 
that 



(HI) 



close to a QNM one can try to iterate for the zero of A^ n using a standard scheme such as Muller's method. This 
strategy works fine for rapidly damped modes (and typically also for the fluid f-mode), but does usually not provide a 
reliable estimate for an imaginary part that is several orders of magnitude smaller than the real part. An alternative 
method is inspired by resonant scattering problems in quantum theory. Letting ui n = a n + i/3 n (with a n and /3 n both 
positive) we can easily identify the real part a n from the position of the minimum of the standard "Breit-Wigner 
resonance" in a graph of log |Ai n |, cf. the example shown in Fig. ^. It is not all that simple to extract the imaginary 
part, however. To do this one must approximate the half- width of the peak, eg. by curve fitting. To achieve satisfactory 
precision in this process is not trivial. 
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FIG. 6: A graph of the incoming amplitude Ai n vs the real part of the frequency. In the left panel we see a standard "Breit- 
Wigner resonance." The right panel illustrates the properties of the phase of the ration of the asymptotic amplitudes. We use 
the location of the two poles and zeros to deduce the frequency and damping rate of a long-lived pulsation mode. 



To devise a better scheme for determining the damping rate of a very long-lived QNM, we employ the two general 
properties we deduced earlier. From these it is clear that, if we have a zero of A{ n close to the real frequency axis, 
there must also be a zero of ^4 ou t on the opposite side of the axis. How does that alter the above results? Close to a 
zero of A out we will have 



A out (u) ~ (u)-u) n ) 



(112) 



Then we immediately see that the absolute value of the ratio between the asymptotic amplitudes, k, is a smooth, 
slowly varying function of the frequency. But if we consider its phase, we find some interesting and useful features. 
We get (for real uj) 



cj - a n + i(3 n 

0n 



uj — a n 

where 7 = j r + fy, is a complex "constant," and therefore 

Im K _ jj[(uj - a n ) 2 - Pi] + 2^ r f3 n {uj - a n ) 
Re k j r [(ui - a n ) 2 - PI] - 2jil3 n (uj - a n ) 



(113) 



(114) 
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From this we see that, instead of having a singularity at u> = a n we now have a function with two zeros and two 
poles on the real frequency axis, cf. Fig. g. Provided that the calculation of the asymptotic amplitudes can be done 
with sufficiently high precision, the location of these zeros (2:1,2) and poles (^1,2) can readily be deduced. Given this 
information one can show that 

On PS ■ (115) 

P1+P2-Z1- Z 2 

Pn ~ K(zi + z 2 ) - o? n - ziz 2 ] 1/2 (116) 

T r ~ ^Pn 1 ro i / 117 \ 

— ~ — o — = ^ 2 "n - zi - z 2 ] . (117) 

It should be pointed out that one gets four equations for the real and imaginary parts of the QNM frequency, as well 
as the ratio "fr/li- The last equality above can therefore be used as a "sanity check" on the calculation. In essence, 
it provides information about the accuracy of the obtained quantitites. 

In practise, this method works very well even when the imaginary part of the frequency is more than eight orders 



of magnitude smaller than the real part, cf. the results in Table [II. 
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